function [mn,stdv,cc,freqbins,midbins]=stats(alpha,NBins)

global NP Z NZ NDRAWS CrossCorr COEF NV BETAS interaction discountfun

%w=createw(alpha); %NPxNDRAWS
w=sum(bsxfun(@times,Z,reshape(alpha,[1,1,NZ])),3);  %NPxNDRAWS
w(w>500)=500;  %As a precaution to prevent machine infty when exponentiate
w=squeeze(exp(w));
w=w./repmat(sum(w,2),1,NDRAWS); %NPxNDRAWS

mn=zeros(NV,1);
v=zeros(NV,1);
cc=zeros(NV,NV);
freqbins=zeros(NV,NBins);
midbins=zeros(NV,NBins);
wtsum=sum(sum(w,1),2);
ww=reshape(w,NP*NDRAWS,1);
for r=1:NV
    thiscoef=squeeze(BETAS(:,r,:)); %NPxNDRAWS
    mn(r,1)=sum(sum(thiscoef.*w,1),2)./wtsum;
    v(r,1)=sum(sum(((thiscoef-mn(r,1)).^2).*w,1),2)./wtsum;
    cc(r,r)=v(r,1);
    if r>2 & CrossCorr==1 & interaction==0 & discountfun~="QuasiHyper";
      for thisc=3:NV;
         thiscoef2=squeeze(BETAS(:,thisc,:)); %NPxNDRAWS
         thismn2=sum(sum(thiscoef2.*w,1),2)./wtsum;
         cc(r,thisc)=sum(sum(((thiscoef-mn(r,1)).*(thiscoef2-thismn2)).*w,1),2)./wtsum;
         cc(thisc,r)=cc(r,thisc);
      end
    elseif r>3 & CrossCorr==1 & interaction==0 & discountfun=="QuasiHyper";
      for thisc=4:NV;
         thiscoef2=squeeze(BETAS(:,thisc,:)); %NPxNDRAWS
         thismn2=sum(sum(thiscoef2.*w,1),2)./wtsum;
         cc(r,thisc)=sum(sum(((thiscoef-mn(r,1)).*(thiscoef2-thismn2)).*w,1),2)./wtsum;
         cc(thisc,r)=cc(r,thisc);
      end
    elseif r>5 & CrossCorr==1 & interaction~=0 & discountfun=="QuasiHyper";
      for thisc=6:12;
         thiscoef2=squeeze(BETAS(:,thisc,:)); %NPxNDRAWS
         thismn2=sum(sum(thiscoef2.*w,1),2)./wtsum;
         cc(r,thisc)=sum(sum(((thiscoef-mn(r,1)).*(thiscoef2-thismn2)).*w,1),2)./wtsum;
         cc(thisc,r)=cc(r,thisc);
      end
    end
    [histw, midpt] = histwgt(reshape(thiscoef,NP*NDRAWS,1), ww, COEF(r,1), COEF(r,2),NBins); 
    freqbins(r,:)=histw./sum(histw);
    midbins(r,:)=midpt;
end
stdv=sqrt(v);